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Abstract 



This article presents an algorithm that allows modeling of biological networks in a qualitative frame- 
work with continuous time. Mathematical modeling is used as a systems biology tool to answer biological 
questions, and more precisely, to validate a network that describes biological observations and to predict 
Q the effect of perturbations. 

There exist two major types of mathematical modeling approaches: (1) Quantitative modeling, de- 
scribing various chemical species concentrations by real numbers, mainly based on differential equations 
and chemical kinetics formalism; (2) Qualitative modeling, describing chemical species concentrations or 
activities by a finite set of values, the most used being the Boolean modeling. 

One major critics of qualitative modeling is that time in this approach is represented by discrete steps. 
Qualitative approaches, on the other hand, can efficiently describe and predict stable states, but remain 
inconvenient in describing transient kinetics leading to these states. In order to handle these transient 
events, several attempts were made to introduce continuous time in qualitative modeling. They consist in 
refining the time evolution in this formalism by introducing priority classes (some variables are updated 
with a higher priority), variable time delays, etc. 
■ * Here, we propose a modeling approach that is intrinsically continuous in time. The algorithm pre- 

sented here fills the gap between qualitative and quantitative modeling. It is based on continuous time 
Markov process applied on a Boolean state space. In order to describe the temporal evolution, we ex- 
plicitly specify the transition rates for each node. For that purpose, we built a language that can be seen 
as a generalization of Boolean equations. The values of transition rates have a natural interpretation: it 
. is the inverse of the time for the transition to occur. Mathematically, this approach can be translated 

in a set of ordinary differential equations on probability distributions; therefore, it can be seen as an 
approach in between quantitative and qualitative. 

We developed a CH — h software, MaBoSS, that is able to simulate such a system by applying Ki- 
netic Monte-Carlo (or Gillespie algorithm) in the Boolean state space. This software, parallelized and 
optimized, computes temporal evolution of probability distributions and can also estimate stationary 
distributions. Applications of Boolean Kinetic Monte-Carlo have been demonstrated for two qualita- 
tive models: a toy model and a published p53/Mdm2 model. Our approach allows to describe kinetic 
phenomena which were difficult to handle in the original models. In particular, transient effects are 
represented by time dependent probability distributions, interpretable in terms of cell populations. 



Background 

Mathematical models of signaling pathways can be seen as tools to answer biological questions. The most 
widely used mathematical formalisms to answer these questions are ordinary differential equations (ODEs) 
and Boolean modeling. 
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Ordinary differential equations (ODEs) have been widely used to model signaling pathways. It is the most 
natural formalism for translating detailed reaction networks into a mathematical model. Indeed, equations 
can be directly derived using mass action laws, Michaelis-Menten kinetics or Hill functions for each reaction 
in order to account for the observed behaviors. This framework has limitations, though. The first one 
concerns the difficulty of assigning the kinetic parameter values in the model. Ideally, these parameters 
would be extracted from experimental data. However, they are often chosen so as to fit qualitatively the 
expected phenotypes. The second limitation arises when studying cell population heterogeneity. In this 
case, ODEs are no longer appropriate since the approach is deterministic and thus focuses on the average 
behavior. To include non-determinism, an ODE model needs to be transformed into stochastic chemical 
model. In this formalism, a master equation is written on the probabilities of number of molecules for each 
species (instead of being written on ODEs of continuous concentrations). In the translation process, the 
same parameters used in ODEs (more particularly in ODEs written with mass action law) can be used in 
the master equation, but in this case, the number of initial conditions explodes along with the computation 
time. 

Boolean formalism is another formalism used to model signaling pathways where genes/proteins are 
parameterized by Os and Is only. It is the most natural formalism to translate an influence network into a 
mathematical model. In such networks, each node corresponds to a species and each arrow to an interaction 
or an influence (positive or negative). In a Boolean model, a logical rule integrating the signs of the input 
links is assigned to each node. As a result, there are no real parameter values to adjust besides choosing the 
appropriate logical rules that best describe the system. In this paper, we will refer to a network state as a 
state in which each node of the influence network has a Boolean value. The set of all possible transitions 
between the network states is defined as a transition graph. There are two types of transition graphs, one 
deduced from the synchronous update strategy [l] , for which all the nodes that can be updated are updated 
in one transition, and another one deduced from the asynchronous update strategy [2], for which only one 
node is updated in one transition. In the Boolean formalism, each transition can be interpreted as a time 
step. Real characterization of biological time is not taken into consideration. However, in many biological 
problems, time plays a crucial role. 

As mentioned for ODE, stochasticity is an important aspect when studying cell populations, It can be 
done: on nodes (by randomly flipping the node state [3|[4] ) , on the logical rules (by allowing to change an 
AND gate by an OR gate [5]), and on the update rules (by defining the probability and the priority of 
changing one Boolean value in an asynchronous strategy [6] or by adding noise to the whole system in a 
synchronous strategy [7]). One of the main drawbacks of the Boolean approach is the explosion of solutions. 
In an asynchronous update strategy, the transition graph can reach 2# nodes . 

Both discrete and continuous frameworks have advantages and disadvantages above-mentioned. We pro- 
pose here to combine some of the advantages of both approaches with what we call the "Boolean Kinetic 
Monte-Carlo" algorithm (BKMC). It consists of a natural generalization of the asynchronous Boolean dy- 
namics [2], with a direct probabilistic interpretation. In BKMC framework, the dynamics is parameterized 
by a biological time, the order of update is noisy (less strict than priority classes introduced in GINsim [8]). A 
BKMC model is specified by logical rules as in regular Boolean models but with a more precise information: 
a numerical rate is added for each transition. 

BKMC is best suited to model signaling pathways in the following cases: 

• The model is based on an influence network, because BKMC is a generalization of the asynchronous 
Boolean dynamics. Typical examples are published models of cell cycle [6j. 

• The model describes processes for which information about biological time is known (relative rate, 
order, etc.), because in BKMC, time is parameterized by a real number. This is typically the case 
when studying developmental biology, where animal models provide time changes of gene/protein 
activities ffl. 

• The model describes heterogeneous cell population behavior, because BKMC has a probabilistic inter- 
pretation. For example, modeling heterogeneous cell population can help to understand tissu formation 
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based on cell differentiation 
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The model can contain many nodes (up to 64 in the present implementation), because BKMC is a 
simulation algorithm that converges fast. This can be useful for big models that have already been 
modeled with a discrete time Boolean method 11 , in order to obtain a finer description of transient 



effects. 

Previous published works have also introduced a continuous time approach in the Boolean frarnework( 
15 , 16 , 17 , 18 ). We will review them and present BKMC approach. We will describe the 



12 , 13 



14 



C++ software, MaBoSS, developed to implement BKMC algorithm and illustrate its use with two examples, 
a toy model and a published model of p53-MDM2 interaction. 

All abbreviations, definitions and algorithms used in this article can be found in Supplementary Material. 
Throughout the article, all terms that are italicized are defined in the Supplementary Material (Part 3). 



Results and discussion 

BKMC for continuous time Boolean model 
Continuous time in Boolean modeling: past and present 

In our context, we briefly recall that, in Boolean approaches for modeling networks, we define the state of 
each node of the network by a Boolean value (node state) and the network state by the set of node states. 
Any dynamics in the transition graph is represented by sequences of network states. A node state is based 
on the sign of the input arrows and the logic that links them. The dynamics can be deterministic (in the 
case of synchronized update [1 ) or non-deterministic (in the case of asynchronized update [2] or probabilistic 
Boolean networks [7]). 

The difficulty to interpret a dynamics in terms of biological time has led to several works that have 
generalized Boolean approaches. These approaches can be divided in two classes that we call explicit and 
implicit time for discrete steps. 

The explicit time for discrete steps consists of adding a real parameter to each node. These parameters 
correspond to the time associated to each node state before it flips to another one. They need to be set for 
each node state ( [13], [l2])- Because data about these time lengths are difficult to extract from experimental 
studies, some works have included noise in the definition of these time parameters [18] . The drawback of 
this method is that the computation of the Boolean model becomes sensitive to both the type of noise and 
the initial conditions. As a result, these time parameters become new parameters that need to be tuned 
carefully and thus add complexity to the modeling. 

The implicit time for discrete steps consists of adding probabilities to each transition of the transition 
graph, in the case of non-deterministic transitions. It is argued that these probabilities could be interpreted 
as adding time to a biological process. As an illustration, let us assume a small network of two nodes, A and 
B. At time t, A and B are inactive: [AB] = [00]. In the transition graph, there exist two possible transitions 
at t+1: [00] — > [01] and [00] — > [10]. If the first transition has a significant higher probability than the second 
one, then we can conclude that B has a higher tendency to activate before A. Therefore, it is equivalent 
to say that the activation of B is faster than the activation of A. Thus, in this case, the notion of time 
is implicitly modeled by setting probability transitions. In particular, priority rules, in the asynchronous 
strategy, consist of putting some of these probabilities to zero [6 . In our example, if B is faster than A 
then the probability of the transition [00] — >• [10] is zero. As a result, the prioritized nodes always activate 
before the others. From a different perspective but keeping the same idea, Vahedi and colleagues [14] have 
set up a method to deduce explicitly these probabilities from the duration of each discrete step. With the 
implementation of implicit time in a Boolean model, the dynamics remains difficult to interpret in terms of 
biological time. 

As an alternative to these approaches, we propose BKMC algorithm. 
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Properties of BKMC algorithm 



BKMC algorithm was built such as to meet the following principles: 

• The state of each node is given by a Boolean number (0 or 1) (referred to as node state); 

• The state of the network is given by the set of node states (referred to as network state); 

• The update of a node state is based on the signs linking the incoming arrows of this node and the 
logic; 

• Time is represented by a real number; 

• Evolution is stochastic. 

For that, we choose to describe the time evolution of network states by a Markov process with con- 
tinuous time. Therefore, the dynamics is defined by transition rates inserted in a master equation (see 
Supplementary Material, section 1.1). Transitions for which a rate is specified are based on asynchronous 
Boolean dynamics. 



Markov process for Boolean model 

In BKMC, we adapt the Markov process to the Boolean approach. 

Consider a network of n nodes (or agents, that can represent any species, i.e. mRNA, proteins, complexes, 
etc.). In a Boolean framework, the network state of the system is described by a vector S of Boolean values, 
i.e. Si € {0, 1}, i = 1, . . . , n where Si is the state of the node i. The set of all possible network states, also 
referred to as the network state space, will be called S. 

A stochastic description of state evolution is represented by a stochastic process s(t) defined on t € I C R 
to the network state space, where I is an interval: for each time t G I C M, s(t) represents a random variable 
to the network state space: 

P [s(t) = S] € [0,1] for any state SeE 
£P[*(f) = S] = 1 (1) 

S6S 

Notice that the random variables s(t) are not independent, therefore P [s(t) = S, s(i') = S'] ^ P [s(t) = S] P [s(tf) = S'] 
From now on, we define P [s(t) = S] as instantaneous probabilities (or first order probabilities). Since the 
instantaneous probabilities do not define the full stochastic process, all possible joint probabilities should 
also be defined. 

In order to simplify the stochastic process, Markov property is imposed. It can be expressed in the 
following way: "the conditional probabilities in the future, related to the present and the past, depend only 
on the present" (see Supplementary Material, section 1.1 for the mathematical definition). A stochastic 
process with the Markov property is called a Markov process. 



Any Markov process can be defined by (see Van Kampen 19 , chapter IV): 

1. an initial condition: 

P[.s(0) = S] ;VS G £ (2) 

2. conditional probabilities (of a single condition): 

P [s(t) = S\s(t') = S'] ; VS, S'eE; W, t e I;t' <t (3) 



Concerning time, two cases can be considered: 
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Time is discrete: t £ I = {to,ti, ■ ■ •}. In that case, it can be shown [20] that all possible conditional 
probabilities are function of transition probabilities: P [s(tj) = S|s(ij_i) = S']. In that case, a Markov 
process is often named a Markov chain. 



Time is continuous: t £ I = [a, b]. In that case, it can be shown 19 that all possible conditional 
probabilities are function of transition rates: P(s'->-s)W <= [0)°°[ ( sce Supplementary Material, section 
1.1 for the definition of transition rates). 

Notice that a discrete time Markov process can be derived from continuous time Markov process, called 
Jump Process, with the following transition probabilities: 

" S^S' 



If the transition probabilities or transition rates are time independent, the Markov process is called a 
time independent Markov process. In BKMC, only this case will be considered. For a time independent 
Markov process, the transition graph (often called Boolean state graph in the Boolean framework) can be 
defined as follows: a transition graph is a graph in S, with an edge between S and S' if and only if ps-s-s' > 
(or P [s(U) = S\s(U-i) = S'] > if time is discrete). 

Asynchronous Boolean dynamics as a discrete time Markov process 

Asynchronous Boolean dynamics (2 is widely used in Boolean modeling. It can be easily interpreted as 



discrete time Markov process 21.22 as shown below. 

In the case of asynchronous Boolean dynamics, the system is given by n nodes (or agents), with a set 
of directed arrows linking these nodes and defining a network. For each node i, a Boolean logic B*(S) is 
specified that depends only on the nodes j for which there exists an arrow from node j to i (e.g. B\ = 
S3 AND NOT 54, where £3 and £4 are the Boolean values of nodes 3 and 4 respectively). The notion of 
asynchronous transition (AT) can be defined as a pair of network states (S, S') £ E, written (S — > S') such 
that 

Sj = Bj(S) for a given j 

S[ = Siiovi^j (4) 

To define a Markov process, the transition probabilities P [s{tj) = S|s(tj_i) = S'] can be defined such 
that: given two network states S and S', let 7(S) be the number of asynchronous transitions from S to all 
possible states S'. Then 

P [s{ti) = S'|s(^_i) = S] = 1/t(S) if (S -> S') is an AT 

P [s(U) = S'|s(ii_i) = S] = if (S -> S') is not an AT (5) 

In this formalism, asynchronous Boolean dynamics completely defines a discrete time Markov process 
when the initial condition is specified. Notice that here the transition probabilities are time independent, 
i.e. P [s(U) = S\s(ti-i) = S'] = P [s(t l+1 ) = S\s(U) = S'}. Other definition of 7(S) can be used, defining 
another Markov process that have the same transition graph. Therefore, the approaches mention above, 
that introduce time implicitly by adding probabilities to each transition of the transition graph, can be seen 
as a generalization of the definition of 7(S). 

Continuous time Markov process as a generalization of asynchronous Boolean Dynamics 

To transform the discrete time Markov process described above in a continuous time Markov process, tran- 
sition probabilities should be replaced by transition rates P(s^s')- In that case, conditional probabilities are 
computed by solving a master equation (equation 2 in Supplementary Material, section 1.1). We present 
below the corresponding numerical algorithm (Kinetic Monte-Carlo |23|). 
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Because we want a generalization of asynchronous Boolean dynamics, transition rates P(s-j-s') ar e non- 
zero only if S and S' differ by only one node. In that case, each Boolean logic -Bi(S) is replaced by two 

functions i?" p / down (S) g [0, oo[. The transition rates arc defined as follow: if i is the node that differs from 
S and S', 

P(s-so = iC(S)if^ = 

P(s-S') = i? down (S) ifS 4 = l (6) 

i?" p corresponds to the activation rate of node i, i? down corresponds to the inhibition rate of node 
i. Therefore, the continuous Markov process is completely defined by all these i? u P/ down and an initial 
condition. 



Asymptotic behavior of continuous time Markov process 



In the case of continuous time Markov process, instantaneous probabilities always converge to a station- 
ary distribution (sec Supplementary Material, corollary 2, section 1.2). A stationary distribution of a 
given Markov process corresponds to the set of instantaneous probabilities of a stationary Markov pro- 
cess which has the same transition probabilities (or transition rates) of the given discrete (or continuous) 
time Markov process. A stationary Markov process has the following property: for every joint probability 
P [s(h) = SW,s(i 2 ) = S( 2 ),...] and Vr, 



S (ti) = S«, s (t 2 ) = S« ... 



s(h+r) = S« .s(t 2 + r) = S«, 



(7) 



Notice that instantaneous probabilities P [s(t) = S] of a stationary stochastic process are time indepen- 
dent. 

The asymptotic behavior of a continuous time Markov process can be detailed by using the concept of 
indecomposable stationary distributions: indecomposable stationary distributions are stationary distribu- 
tions that cannot be expressed as linear combination of different stationary distributions. Notice that a 
linear combination of stationary distributions is also a stationary distribution, up to a constant. This comes 
from the fact that instantaneous probabilities are solutions of a master equation, which is linear (see Sup- 
plementary Material, equation 2, section 1.1). Therefore, a complete description of the asymptotic behavior 
is given by the linear combination of indecomposable stationary distributions, to which the Markov process 
converges. 



Oscillations and cycles 

In order to describe a periodic behavior, the notion of cycle and oscillation for a continuous time Markov 
process is defined precisely. 

A cycle is a loop in the transition graph. This is a topological characterization that does not depend on 
the exact value of transition rates. It can be shown that a cycle with no outcoming edge corresponds to an 
indecomposable stationary distribution (see Supplementary Material, corollary 1, section 1.2). 

The question is then to link the notion of cycle to that of periodic behavior of instantaneous probabilities. 
Instantaneous probabilities cannot be perfectly periodic; at most, they have a damped ocscillating behavior 
(see Supplementary Material, section 1.3). Let us define formally a damped oscillatory Markov process as a 
continuous time process that has at least one instantaneous probability with an infinite number of extrema. 

According to theorems described in Supplementary Material (theorems 6-8 and Corollary 3, section 
1.3), a necessary condition for having damped oscillation is that the transition matrix (see Supplementary 
Material, equation 4, section 1.1) has at least one non-real eigenvalue. In that case, there always exists an 
initial condition that produces damped oscillations. For the transition matrix to have a non-real eigenvalue, 
a Markov process needs to have a cycle. However, the reverse is not true. In the toy model of single cycle, 
presented in the section of examples, non-real eigenvalues may or may not exist, according to different sets 
of transition rates, although the transition graph remains the same. 
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BKMC: Kinetic Monte-Carlo (Gillespie algorithm) applied to continuous time asynchronous 
Boolean Dynamics 

It has been previously stated that a continuous time Markov process is completely defined by its initial 
condition and its transition rates. For computing any conditional probability (and any joint probability), a 
set of linear differential equations has to be solved (the master equation). Theoretically, the master equation 
can be solved exactly by computing the exponential of the transition matrix (see Supplementary Material, 
equation 5, section 1.1). However, because the size of this transition matrix is 2" x 2™, practical computation 
soon becomes impossible if n is large. To remedy this problem, it is possible to use a simulation algorithm 
that samples the probability space by computing time trajectories in the network state space. 

The Kinetic Monte-Carlo [23] (or Gillespie algorithm [24]) is a simple algorithm for exploring the prob- 
ability space of a Markov process defined by a set of transition rates. In fact, it can be understood as a 
formal definition of a continuous time Markov process. This algorithm produces a set of realizations or 
stochastic trajectories of the Markov process, given a set of uniform random numbers in [0, 1[. By definition, 
a trajectory S(t) is a function from a time window [0,i max ] to S. The set of realizations or stochastic tra- 
jectories represents the given Markov process in the sense that these trajectories can be used to computed 
probabilities. Practically, a finite set of these trajectories is produced, then probabilities are estimated from 
this finite set (as described below). The algorithm is based on an iterative step: from a state S at time to 
(given two uniform random numbers), it produces a transition time St and a new state S', with the following 
interpretation: the trajectory S(t) is such that S(t) = S for t E [to, to + St[ and S(f + St) = S'. Iteration of 
this step is done until a specified maximum time is reached. The initial state of each trajectory is based on 
the (probabilistic) initial condition, that needs also to be specified. 

The exact iterative step is the following, given S and two uniform random number u, v! € [0, 1[: 

1. Compute the total rate of possible transitions for leaving S state: 

Ptot = J2s' P(S-J-S')- 

2. Compute the time of the transition: St = — log(tt)/ptot 

3. Order the possible new states S'^\j = 1 . . . and their respective transition rates p& = P(S->-S'0))- 

4. Compute the new state S'( fc ) such that X^=o Pi < W Ptot) < Y^j=o Pi (°y convention, p^ = 0). 

The application of this algorithm to continuous time Markov process in network state space will be 
referred to as Boolean Kinetic Monte-Carlo or BKMC. 

Practical use of BKMC, through MaBoSS tool 

Biological data are summarized into an influence network with logical rules associated to each node of 
the network. The value of one node depends on the value of the input nodes. For BKMC, another layer 
of information is provided when compared to the standard definition of Boolean models: transition rates 
are provided for all nodes, specifying the rates at which the node turns on and off based on their logic 
for both the on and off rules. This refinement conserves the simplicity of Boolean description but allows 
to reproduce the observed biological dynamics. The parameters do not need to be exact as in nonlinear 
ordinary differential equation models but they can be used to illustrate the relative speed of reactions. For 
that purpose, we developed a software tool, MaBoSS, that applies BKMC algorithm. MaBoSS stands for 
Markov Boolean Stochastic Simulator. Practically, MaBoSS needs two input files: one describing the network 
and its transition rates, one describing the parameters controlling the different estimates described in the 
Methods section. Source code, reference card and examples are available on the web: https://maboss.curie.fr. 

Transition rates with MaBoSS language: modeling biological processes 

MaBoSS defines transition rates P(s^s') by the functions ^P/ down (g) (see equation pi. The format of these 
functions is very flexible. It includes all Boolean operators (AND, OR, NOT, XOR), arithmetic operators (+,- 
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,* /), external variables, node variables, comparison operators and the conditional operator (?:). Examples 
of the use of the language are given below to illustrate three different cases: different speeds for different 
inputs, buffering effect and the translation of discrete variables (with more than the 2 values, and f ) in 
MaBoSS. 

• Modeling different speeds of incoming influences: suppose that C is activated by A and B, but that B 
can activate C faster than A. In this case, we write: 



node C { 

rate_up= B ? $kb : (A ? $ka : 0.0); 
rate_down= (A I B ) ? . : 1 . 0} 



When C is off (equal to 0), C is activated by B at a speed kb. If B is absent, then C is activated by A 
at a speed Ska. If both are absent, C is not activated. Note that if both A and B are present, because 
of the way the logic is written in this particular case, C is activated at a speed $kb. When C is on 
(equal to 1), C is inactivated at a rate equal to I if A and B are both absent. 

To implement the synergetic effect of A and B, i.e. when both A and B are on, C activates at a rate 
$kab, then we can write: 

node C { 

rate_up= (A & !B ? $ka : 0)+(B k !A ? $kb : 0) + (A & B ? $kab : 0.0); 
rate_down= (A I B ) ? . : 1 . 0} 



• Modeling buffering effect: suppose that B is activated by A, but that B can remain active a long time 
after A has shut down. For that, it is enough to define different speeds of activation and inhibition: 



node B { 

rate_up= A ? 2.0 : 0.0; 
rate_down= A ? 0.0 : 0.001;} 



B is activated by A at a rate equal to 2. When A is turned off, B is inactivated more slowly at a rate 
equal to 0.001. 

• Modeling more than two discrete states for a given node: Suppose that B is activated by A, but if the 
activity of A is maintained, B can reach a second level. For this, we define a second node B Ji (for "B 
high") with the following rules: 



node B { 

rate_up= A ? 1.0 : 0.0; 
rate_down= (A I B_h) ? 0.0 : 1.0;} 

node B_h { 

rate_up= (A & B) ? 1.0 : 0.0; 
rate_down= (A) ? 0.0 : 1.0;} 



In this example, B is separated in two variables: B which corresponds to the first level of B and B_h 
which corresponds to the higher level of B. B is activated by A at a rate equal to 1. If A disappears 
before B has reached its second level B_h then B is turned off at a rate equal to 1. If A is maintained 
and B is active, then BJi is activated at a rate equal to 1. When A is turned off, B_h is inactivated at 
a rate equal to 1. 
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Simulation input parameters in MaBoSS 



To simulate a process in MaBoSS, a set of parameters need to be adjusted (see simulation parameters 
in reference card). MaBoSS assigns default values, however, they need to be tuned for each model to 
achieve optimal performances: best balance between the convergence of estimates and the computation 
time. Therefore, several simulations should be run with different set of parameters for best tuning. 

• Internal nodes: node. is -internal As explained in Methods ("Initial conditions and outputs"), internal 
nodes correspond to species that are not measured explicitly. Practically, the higher the number of 
internal nodes, the better the convergence of the BKMC algorithm. 

• Time window for probabilities: timetick 

This parameter is used to compute estimates of network state probabilities (see "Network state prob- 
abilities on time window" in Methods). A time window can be set as the minimum time needed for 
nodes to change their states. This parameter also controls the convergence of probability estimates. 
The larger the time window parameter, the better the convergence of probability estimates. With 
practice, the tradeoff between timetick parameter value and the convergence speed will be defined. 

• Maximum time: maxJime 

MaBoSS produces trajectories for a predefined amount of time set by the parameter: max.time. This 
maximum time needs to be specified. If the time of the biological process is known, then the maximum 
time parameter can be explicitly set. If the time of the biological process is not known, then there exists 
a more empirical way to set the maximum time. It is advised to choose a maximum time parameter 
that is slightly bigger than the inverse of the smallest transition rate. 

Note that the computing time in MaBoSS is proportional to this maximum time. Moreover, the choice 
of the maximum time impacts the stationary distribution estimates (see below): a longer maximum 
time increases the quality of these estimates. 

• Number of trajectories: sample-count 

This parameter directly controls the quality of BKMC estimation algorithm. Practically, the conver- 
gence of the estimates increases as the number of trajectories is increased. 

• Number of trajectories (statdist-traj-count) and similarity threshold (statsdisLcluster-threshold) for 
stationary distribution estimates 

The statdist-traj-count parameter corresponds to a subset of trajectories use only for stationary dis- 
tribution estimates (the statdist-traj -count first trajectories are chosen by the algorithm). To avoid 
explosion of computing time, this parameter needs to be lower than the number of trajectories (rather 
than equal to). 

The statsdisLcluster -threshold parameter corresponds to the threshold for constructing the clusters of 
stationary distribution estimates. Ideally, it should be set to a high value (close to 1). However, if the 
threshold is too high then the clustering algorithm might not be efficient. 

For optimal results, the identification of the full set of indecomposable stationary distributions should 
be done in the following way: 

— Run a simulation with an initial condition and maximum time, with a reasonable similarity thresh- 
old (around 0.8) and a number of trajectories (around 1000). As a response, MaBoSS provides 
a set of indecomposable stationary distributions (it corresponds to the stationary distributions 
associated to each cluster). 

— Select the states with non-zero probability in the set of indecomposable stationary distributions 
and set them as initial conditions, increase the maximum time (max-time), the number of trajec- 
tories {statdist-traj -count) and/or the similarity threshold (statsdisLcluster-threshold) . 
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— Run the simulations for these new parameters. 

— If new indecomposable stationary distributions appear, start again the two previous steps. Stop 
when the indecomposable stationary distributions remain stable with respect to the simulation 
parameters, i.e. after several rounds. 

Comparison with biological data 

The relationship between a model and experimental data is strongly dependent on the type of model and 
the question that needs to be solved. 

Because MaBoSS is based on Boolean modeling, the biological data need to be discretized. Each node of 
the model should represent discrete levels of the respective species (mRNA, protein, protein complex, etc.). 
It is possible to have more than two discrete levels in a model, as shown in the example "Modeling more 
than two discrete states for a given node" . 

The transitions rates are positive numbers that should be introduced in a model; it is possible to extract 
them from experimental data, using the following property: the rate of a given transition is the inverse of 
the mean time, for this transition to happen. It should be noticed than BKMC is an algorithm based on a 
linear equation (Supplementary Material, equation 2, section 1.1); therefore, small variations of transition 
rates won't affect qualitative behavior of a model. 

The basic outputs of BKMC algorithm are network state probabilities over time. These can be interpreted 
in terms of a cell population, once experimental data are discretized. The asymptotic behavior of a model, 
represented by a linear combination of indecomposable stationary distributions, can be interpreted as a 
combination of cell sub-populations. More precisely, consider an indecomposable stationary distribution and 
the associated set of network states with non-zero probability. A cell in such a network state can only evolve 
in other network states with non-zero probability, within the same indecomposable stationary distribution 
(Supplementary Material, corollary 1, section 1.2 and by using the definition of strongly connected component 
with no outcoming edge). Therefore, the set of network states with non-zero probability can be interpreted 
as a sub-population whose cells evolve only within the sub-population. 

Examples 

Two models using BKMC applied to Boolean networks are given as examples. The first one is a toy model, 
illustrating the dynamics of a single cycle. The second one uses a published Boolean model of p53-Mdm2 
response to DNA damage. Note that MaBoSS has been used for these two examples, but Markov process 
can be computed directly, without our BKMC algorithm because the model is small enough (by computing 
exponential of transition matrix, see Supplementary Material, section 1.1), as proposed in [16]. BKMC is 
best suited for larger networks, when the network state space is too large to be computed with standard 
existing tools (>~ 2 10 ). These examples were chosen for their simplicity, and because they illustrate how 
global characterizations (entropy and transition entropy, see "Entropies" in Methods) can be used. 

For the purpose of this article, we built the transition graphs for both examples (with GINsim [8j) in 
order to help the reasoning. However, BKMC algorithm does not construct the transition graph explicitly. 

Toy model of a single cycle 

We consider three species, A, B and C, where A is activated by C and inhibited by B, B is activated by A 
and C is activated by A or B (Figure [T]). 

The model is defined within the language of MaBoSS (defined in the web page, https://maboss.curie.fr). 
The associated transition graph is shown in figure [2j In this model, the only stationary distribution is the 
fixed point [ABC] = [000]. We studied two cases: when the rate of the transition from state [001] to state 
[000] is either fast or slow (inactivation of C). We will refer to this transition rate as the escape rate. 

In the first case, when the escape rate is fast, we set the parameter for the transition to a high value 
(rate_up = 10). In figure |3j we notice that the probability that [ABC] is equal to [000] converges to 1. We 
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// Toy model of a 4 states cyde 


$Au=l; 


NodeC 


$Bu=2; 


{ 


$Bd=B; 


rate up=0.0; 


$Ad=4; 


rate dowrH(NOTA) AND (NOT B)) ? $escape : 0.0 


$escape=10; // self-degradation of C that makes the 4cyde unstable. 




Cistate=l; 


'} 


Aistate=l; 




B.istate=l; 


Node A 


Arefstate=0; 


{ 


discretetime -0; 


rate_up={TAND(NOTB)) ? $Au : 0.0; 


use physrandgen =FAL£E; 


rate down=B ? $Ad : 0.0; 


seedpseudorandom -100; 


} 


samplecount =500000; 




max_time=5; 


NodeB 


tinre_tid<=0.01; 


{ 


threadcount =4; 


rate_up=A?$Au;0.0; 


statdist traj count =100; 


rate down =A? 0.0 :$Ad; 


statdist duster threshold =0.9; 


} 





Figure 1: Toy model of a single cycle. 
(C) Simulation parameters 



(A) Influence network. (B) Logic and transition rates of the model. 




Figure 2: Transition graph for the toy model. The node state should be read as [ABC] = [***]. [ABC] = [100] 
corresponds to a state in which only A is active. The nodes in green belong to a cycle, the node in red is 
the fixed point and the other state nodes are in blue. 



can conclude that [ABC] = [000] is a fixed point. In addition, the entropy and the transition entropy converge 
to zero. With BKMC, these properties correspond to a signature of a fixed point. The peak in the entropy 
(between times and 0.6) corresponds to a set of transiently activated states before reaching the fixed point. 

In the second case, when the escape rate is slow, we set the parameter for the transition to a low value 
(rate_down = 10~ 5 ). The model seems to show another stationary distribution. As illustrated in figure [ij 
the transition entropy is and remains close to zero but the entropy does not converge to zero, which is the 
signature of a cyclic stationary distribution (see "Entropies" in Methods). This corresponds to the cycle 
[111] — > [Oil] — > [001] — > [101] in the transition graph However, as seen in the transition graph, one 
state in the cycle will eventually lead the fixed point (through the transition [001] — > [000] in figure |2|. 
Therefore, if the temporal evolution is plotted on a larger time scale (Figure [5]), it looks similar to the case 
of fast escape rate. This case can be anticipated. Indeed, the value of the transition entropy of figure [4] is 
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Figure 3: BKMC algorithm applied to the toy model, with a fast escape rate. Time trajectory of probabilities 
([ABC] = [000] and [ABC] = [1**] where * can be either or 1), the entropy (H) and the transition entropy 
(TH) are plotted. Because the probability of [ABC] = [000] converges to 1, [ABC]=[000] is a fixed point. The 
asymptotic behavior of both the entropy and the transition entropy is also the signature of a fixed point. 

not exactly zero, but 10~ 4 . Therefore, the cyclic behavior is not stable. We can conclude on stable cyclic 
behaviors only when the transition entropy is exactly zero. 

By considering the spectrum of the transition matrix (see Supplementary Material, section 1.1 and proof 
of theorem 4), it can be proven that the model with a slow escape rate is a damped oscillatory process 
whereas the model with a large escape rate is not a damped oscillatory process. As mentioned previously, 
a cycle in the transition graph may or may not lead to an oscillatory behavior. Moreover, if the transition 
entropy seems to converge to a small value on a small time scale, and the entropy does not, this behavior 
illustrates the case of a transient cycle in the transition graph. 

p53-Mdm2 signaling 

We consider a model of p53 response to DNA damage [l8|. p53 interacts with Mdm2, which appears in two 
forms, cytoplasmic and nuclear. On one hand, p53 upregulates the level of cytoplasmic Mdm2 which is then 
transported into the nucleus and inhibits the export of nuclear Mdm2. On the other hand, Mdm2 facilitates 
the degradation of p53 through ubiquitination. In the model, stress regulates the level of DNA damage, 
which in turn participates in the degradation process of Mdm2. p53 inhibits the DNA damage signal by 
promoting DNA repair. Here, stress is not shown explicitly (Figure [6]). 

The model is defined within the language of MaBoSS, with two levels of p53 as it is done in Abou-Jaoude 
et al. [18] . The model is implemented in MaBoSS (provided in the web page (https://maboss.curie.fr) along 
with the simulation parameters) . The associated transition graph is given in figure [7j It shows the existence 
of two cycles and of a fixed point [p53 Mdm2C Mdm2N Dam] = [0010]. 

In order to represent the activity of p53, time evolution of its expectation value is shown in figure [HJ with 
the initial condition: [p53 Mdm2C Mdm2N Dam] = [0*11]. The activation of p53 seems to be transient. 

The qualitative results obtained with MaBoSS are similar to those of Abou-Jaoude and colleagues. 
However, at the level of cell population, some discrepancies appear: in figure [5J no damped oscillations can 
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Figure 4: BKMC algorithm applied to the toy model, with a slow escape rate. Time trajectory of probabil- 
ities ([ABC] = [000] and [ABC] = [1**]), the entropy (H) and the transition entropy (TH) are plotted. The 
asymptotic behavior of both the entropy and the transition entropy seems to be the signature of a cycle. 



be seen as opposed to figure 8 of their article. The reason is that in their computations, noise imposed on time 
is defined by a square distribution on a limited time frame, whereas in BKMC, Markovian hypotheses imply 
that the noise distribution is more spread out from to infinity. The consequence is that synchronization is 
lost very fast. Damped oscillations could be observed with BKMC with a particular set of parameters: fast 
activation of p53 and slow degradation of p53 (results not shown) . 

With MaBoSS, we clearly interpret the system as a population and not as a single cell. In addition, we 
can simulate different contexts, presented in the initial article as different models, within one model that used 
different parameters to account for these contexts (see web page https://maboss.curie.fr for more details). 

Note that the existence of transient cycles, as shown in the toy model, can be deduced from the time 
trajectory of the entropy that is significantly higher than the time trajectory of the transition entropy (which 
is non zero therefore the transient cycles are not stable) (Figure [9| . An other indirect effect of transient 
cycles is the flat part in p53 standard deviation, in figure [8] 

Conclusions 

In this work, we present a new algorithm, Boolean Kinetic Monte-Carlo or BKMC, applicable to dynamical 
simulation of signaling networks based on continuous time in the Boolean framework. BKMC algorithm is a 
natural generalization of asynchronous Boolean dynamics [2], with time trajectories that can be interpreted 
in terms of biological time. The variables of the Boolean model are biological species and the parameters 
are rates of activation or inhibition of these species. These parameters can be deduced from experimental 
data. We applied this algorithm to two different models: a toy model that illustrates a simple cyclic behavior, 
and a published model of p53 response to DNA damage. 

This algorithm is provided within a freely available software, MaBoSS, that can run BKMC algorithm on 
networks up to 64 nodes, in the present version. The construction of a model uses a specific grammar that 



13 




Figure 5: BKMC algorithm applied to the toy model, with a slow escape rate, plotted on a larger time scale. 
Time trajectory of probabilities ([ABC] = [000] and [ABC] = [1**]), the entropy (H) and the transition entropy 
(TH) are plotted. On a large time scale, the asymptotic behavior of both the entropy and the transition 
entropy is similar to the case of the fast escape rate (figure [3| . 




Figure 6: Boolean model of p53 response to DNA damage 



allows to introduce logical rules and rates of node activation/inhibition in a flexible manner. This software 
provides global and semi-global outputs of the model dynamics that can be interpreted as signatures of the 
dynamical behaviors. These interpretations become particularly useful when the network state space is too 
large to be handled. The convergence of BKMC algorithm can be controlled by tuning different parameters: 
maximum time of the simulation, number of trajectories, length of a time window on which the average of 
probabilities is performed, and the threshold for the definition of stationary distributions clusters. 

The next step is to apply BKMC algorithm with MaBoSS on large signaling networks, e.g. EGFR 
pathway, cell cycle signaling, apoptosis pathway, etc. The translation of existing Boolean models in MaBoSS 
is straightforward but requires the addition of transition rates. In these future works, we expect to illustrate 
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Figure 7: Transition graph of the p53 model. The node states should be read as [p53 Mdm2C Mdm2N Dam] 
_ ^jjgj-g * can ^g gither or 1). For instance, [p53 Mdm2C Mdm2N Dam] = [1000] corresponds to a 

state in which only p53 (at its level 1) is active. The nodes in green and the nodes in light blue belong to 
two cycles, the node in red is the fixed point and the other state nodes are in dark blue. 
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Figure 8: Time trajectories of p53 expectation value and standard deviation 



the advantage of BKMC on other simulation algorithms. Moreover, in future developments of MaBoSS, we 
plan to introduce methods for sensitivity analyses and to refine approximation methods used in BKMC, and 
generalize Markov property. 



Methods 

BKMC generates stochastic trajectories. Here, we describe how we use and interpret them. 
Network state probabilities on time window 

To relate continuous time probabilities to real processes, an observable time window At is defined. A discrete 
time (t € N) stochastic process s(r) (that is not necessary Markovian) can be extracted from the continuous 
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Figure 9: Time trajectories of the entropy (H) and the transition entropy (TH) 



time Markov process: 

1 f (T+l)At 

P [s(t) =S] = —j^ dtV [s(t) = S] (8) 
BKMC is used for estimating P [s(t) = S] as follows: 

1. For each trajectory j, compute the time for which the system is in state S, in the window [rAt, (r 4- 
l)At]. Divide this time by At. Obtain an estimate of P [s(t) = S] for trajectory j, i.e. Pj [s(t) = S]. 

2. Compute the average over j of all Pj [s(t) = S] to obtain P [s(r) = S]. Compute the error of this 
average (^/var(P [s(t) = S])/# trajectories). 

Entropies 

Once P [s(t) = S] is computed, the entropy H(t) can be estimated 

H(r) =-J2 l«g 2 (P [*(r) = S]) P [s(r) = S] (9) 
s 

The entropy measures the disorder of the system. A maximum entropy means that all states have the same 
probability, a zero entropy means that one of the states has a probability of one. The estimation of the entropy 
can be seen as a global characterization (by a single real number) of a full probability distribution. The 
choice of log 2 allows to interpret H(t) in an easier manner: 2 H ^ is an estimate of the number of states that 
have a non-negligible probability in the time window [rAt, (r + l)At]. A more computer-like interpretation 
of H(t) is the number of bits that is necessary for describing states of non-negligible probability. 

The Transition Entropy TH is a finer measure that characterizes the system at the single trajectory level; 
it can be computed in the following way: for each state S, there exists a set of possible transitions S — > S'. 
For each of these transitions, a probability is associated: 

P s ^s< = ? S ^ S ' ■ (10) 

By convention, Ps->S' = if there is no transition from S to any other state. 
Therefore, the transition entropy TH can be associated to each state S: 

TH(S) = - log 2 (P S ^s')Ps^S' (11) 

S' 

Similarly, TH(S) — if there is no transition from S to any other state. The discrete time transition entropy 
TH(t) is defined as: 

TH(t)=J2PHt) = S]TH(S) 
s 
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This transition entropy is estimated in the following way: 

1. For each trajectory j, compute the set ($) of visited states (S) in time window [rAt, (r + l)At] and 
their respective duration (/is)- The estimated transition entropy is 

Tff(r). = ^T#(S)g (12) 

2. Compute the average over j of all TH(t)- to obtain TH(t). Compute the error of that average 

{\Jy&t(TH(t))/# trajectories). 

This transition entropy is a way to measure how deterministic the dynamics is. If the transition entropy is 
always zero, the system can only make a transition to a given state (although time of transitions remains 
stochastic). 

If probability distributions on time window tend to constant values (or tend to a stationary distribution) , 
the entropy and the transition entropy can help characterize this stationary distribution such that: 

• A fixed point has zero entropy and zero transition entropy, 

• A cyclic stationary distribution has non-zero entropy and zero transition entropy. 

Entropy and transition entropy can be considered as "global characterization" of time evolution model: 
for a given time window, they always consist of two real numbers, whatever the size of the network. 

Hamming distance distribution 

The Hamming Distance between two states S and S' is the number of nodes that have different node states 
between S and S': 

HD(S,S') = ^{l-6 Si ,sO (13) 

i 

where 5 is the Kronecker delta (Ssi,S'. = 1 if Sj = S[, Ss it s'- = if 5j ^ S-). Given a reference state S re f, the 
Hamming distance distribution (over time) is given by 

P(HD, t) = 5> [s(t) = S] S HDMD{s . Sraf) (14) 
s 

The estimation of discrete time Hamming distance distribution P(HD, r) is similar to that of stochastic 
probabilities on time window. 

The Hamming distance distribution is a useful characterization when the set of instantaneous probabilities 
is compared to a reference state (S re f). In that case, Hamming distance distribution describes how far this 
set is to this reference state. 

The Hamming distance distribution can be considered as a "semi-global" characterization of time evolu- 
tion: for a given time window, the size of this characterization is the number of nodes (to be compared with 
probabilities on time window, whose size is 2# nodos ). 

Initial conditions and outputs 

Inputs Nodes are defined as the nodes for which the initial condition is fixed. Therefore, each trajectory of 
BKMC starts with fixed values of input nodes and random values of other nodes. 

Internal nodes are nodes that are not considered for computing probability distributions, entropies and 
transition entropies. Output nodes are nodes that are not internal. Technically, probabilities are sum up 
over network states that differ only by the state of internal nodes. These internal nodes are only used for 
generating time trajectories through BKMC algorithm. Usually, nodes are chosen to be internal when the 
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corresponding species is not measured experimentally. Mathematically, it is equivalent to transform the 
original Markov process to a new stochastic process (that is not necessary Markovian) defined to a new 
network state space; this new state space is defined by states of output nodes. This raises the question of 
the transition entropy TH: formally, this notion has only a sense within Markovian processes, i.e. when 
there are no internal nodes. Here, we generalize the notion of transition entropy even in the case of internal 
nodes. Suppose that the system is in state S: 

• If the only possible transitions from state S to any other state consist of flipping an internal node, the 
transition entropy is zero. 

• If there is, at least, one transition from state S to another state that flips a non-internal node, then only 
the non-internal nodes will be considered for computing probabilities in equations |10| In particular, 
Ylai ps-s-S' is computed only on non-internal node flipping events. 

Reference nodes are nodes for which a reference node state is specified and for which the Hamming 
distance is computed. In this framework, a reference state is composed of reference nodes for which the node 
state is known and non-reference nodes for which the node state is unknown. Note that non-reference nodes 
may differ from internal nodes. 



Stationary distribution characterization 



It can be shown (see Supplementary Material, corollary 2, section 1.2) that instantaneous probabilities of 
a continuous time Markov process converge to a stationary distribution. Fixed points and cycles are two 
special cases of stationary distributions. They can be identified by the asymptotic behavior of entropy and 
transition entropy (this works only if no nodes are internal) : 

• if the transition entropy and the entropy converge to zero, then the process converges to a fixed point 

• if the transition entropy and the entropy does not converge to zero, then the process converges to a 
cycle 

More generally, the complete description of the Markov process asymptotic behavior can be expressed as 
a linear combination of the indecomposable stationary distributions. 

A set of finite trajectories, produced by BKMC, can be used to estimate the set of indecomposable 



stationary distributions. Consider a trajectory S(t),t G [0,T],i = 1, • • • ,n. Let Is{t) = ^ss(t) - The 
estimation of the associated indecomposable stationary probability distribution (sq) is done by averaging 
over the whole trajectory: 

1 



P [s = S] 



T 



dth(t) 



(15) 



Therefore, a set of indecomposable stationary distribution estimates can be obtained by a set of trajec- 
tories. These indecomposable stationary distribution estimates should be clustered in groups, where each 
group consists of estimates for the same indecomposable stationary distribution. For that, we use the fact 
that two indecomposable stationary distributions are identical if they have the same support, i.e. the same 
set of network states with non-zero probabilities (shown in Supplementary Material, theorem 2). Therefore, 
it is possible to quantify how similar two indecomposable stationary distribution estimates are. A similarity 



coefficient D(sq\ Sq^) £ [0, 1], given two stationary distribution estimates Sq' and Sq' , is defined 



where 



( 

yS6support(s<, ,) ,s<, j) ) 




£ 



support(sg \ Sq ^) = js such that P 



S'esupport(s{ i ,, ,s^' , ) 



>0 



S' 



(16) 



(17) 
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Clusters can be constructed when a similarity threshold a is provided. A cluster of stationary distributions 
is defined as follows: 

C = {s | 3s' e C s. t. D(s , «o) > (18) 

For each cluster C, an associated distribution sc, that estimates an indecomposable stationary distribu- 
tion, can be defined 

P[ Sc = S] = i^P[, = S] (19) 

Errors on this estimate can be computed by 

Err (P [s c = S]) = A /Vax(P [a = S] , s i C)/|C| (20) 

Notice that this clustering procedure has no sense if the process is not Markovian; therefore, no node will 
be considered as internal. 

List of abbreviations 

BKMC: Boolean Kinetic Monte-Carlo 

AT: Asynchronous transition 

ODEs: Ordinary Differential Equations 

MaBoSS: Markov Boolean Stochastic Simulator 
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